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We perform a systematic study of several models that have been proposed for the purpose of 
understanding the motion of driven interfaces in disordered media. We identify two distinct uni- 
versality classes: (i) One of these, referred to as directed percolation depinning (DPD), can be 
described by a Langevin equation similar to the Kardar-Parisi-Zhang equation, but with quenched 
disorder, (ii) The other, referred to as quenched Edwards- Wilkinson (QEW), can be described by a 
Langevin equation similar to the Edwards- Wilkinson equation but with quenched disorder. We find 
that for the DPD universality class the coefficient A of the nonlinear term diverges at the depinning 
transition, while for the QEW universality class either A = or A — + as the depinning transition is 
approached. The identification of the two universality classes allows us to better understand many 
of the results previously obtained experimentally and numerically. However, we find that some 
results cannot be understood in terms of the exponents obtained for the two universality classes 
at the depinning transition. In order to understand these remaining disagreements, we investigate 
the scaling properties of models in each of the two universality classes above the depinning tran- 
sition. For the DPD universality class, we find for the roughness exponent ap = 0.63 ± 0.03 for 
the pinned phase, and cxm = 0.75 ± 0.05 for the moving phase. For the growth exponent, we find 
f3p = 0.67 ± 0.05 for the pinned phase, and (5m = 0.74 ± 0.06 for the moving phase. Furthermore, 
we find an anomalous scaling of the prefactor of the width on the driving force. A new exponent 
ipM = —0.12 ± 0.06, characterizing the scaling of this prefactor, is shown to relate the values of the 
roughness exponents on both sides of the depinning transition. For the QEW universality class, we 
find that ap ~ a.u = 0.92 ± 0.04 and ftp ~ /3m = 0.86 ± 0.03 are roughly the same for both the 
pinned and moving phases. Moreover, we again find a dependence of the prefactor of the width on 
the driving force. For this universality class, the exponent ipM = 0.44 ± 0.05 is found to relate the 
different values of the local ap and global roughness exponent ac ~ 1.23 ± 0.04 at the depinning 
transition. These results provide us with a more consistent understanding of the scaling properties 
of the two universality classes, both at and above the depinning transition. We compare our results 
with all the relevant experiments. 

PACS numbers: 47.55.Mh 68.35.Fx 



I. INTRODUCTION 

Recently the growth of rough interfaces has witnessed 
an explosion of theoretical, numerical, and experimental 
studies, fueled by the broad interdisciplinary aspects of 
the subject 0-0]. Applications can be so diverse as im- 
bibition in porous media, fluid-fluid displacement, fire 
front motion, and the motion of flux lines in supercon- 
ductors @-|lJ. 

In general, a d-dimensional self-affine interface, de- 
scribed by a single-valued function h(x,t), evolves in a 
(d + l)-dimensional medium. Usually, some form of dis- 
order rj affects the motion of the interface leading to its 
roughening. Two main classes of disorder have been dis- 
cussed in the literature. The first, called thermal or "an- 
nealed," depends only on time. The second, referred to 
as "quenched," is frozen in the medium. Early studies 
focused on time- dependent disorder as being responsible 
for the roughening fl8|-|26|]. Here, we consider in detail 
the effect of quenched disorder on the growth. 

The presence of quenched disorder introduces an inter- 
esting analogy between the motion of driven interfaces in 



disordered media and the theory of critical phenomena. 
The continual motion of the interface requires the ap- 
plication of a driving force F. There exists a critical 
force F c , such that for F < F C) the interface will be- 
come pinned by the disorder after some finite time. For 
F > F c , the interface moves indefinitely with an average 
velocity v(F). This suggests that the motion of driven 
rough interfaces in disordered media can be studied as 
a phase transition — which we shall call the depinning 
transition. The velocity of the interface v plays the role 
of the order parameter, since as F — ► F^~ , v vanishes as 

v~f. (1.1) 

We call 9 the velocity exponent, and / = (F — F c )/F c 
the reduced force (Fig. [I]). 

For F — » F c + , large but finite regions of the interface 
are pinned by the disorder. At the transition, the char- 
acteristic length £ of these pinned regions diverges as 

t,~r v , (i.2) 

where v is the correlation length exponent. 

The added phenomenological richness introduced by 
the quenched disorder leads to an increased difficulty. In 
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fact, the problem of identifying the universality classes 
for interface roughening in the presence of quenched dis- 
order and determining the corresponding sets of critical 
exponents has so far remain unsolved. Experimental and 
numerical measurements of some of the critical exponents 
vary considerably, and the question of the existence of 
universality classes bas been raised J^ |r^ , [27| ~[f7| . The 
problem is made more complex by the disagreement of 
most experimentally measured values with analytical cal- 
culations @-fo|. 

In this paper, we address the question of why the same 
exponents vary in value for different systems. We study 
several numerical models introduced to study interface 
roughening in the presence of quenched disorder, and 
identify two universality classes. We show that one of 
those universality classes can be identified with the case 
studied analytically. The exponents for the other uni- 
versality class cannot be determined from renormaliza- 
tion groups calculations. However, mappings to directed 
percolation (DP) and isotropic percolation allow us to 
estimate most scaling exponents. 

The identification of the universality classes and the 
calculation of the respective scaling exponents at the de- 
pinning transition enable us to understand the results for 
most of the models proposed and for some of the exper- 
iments. We shall see that several experimental results 
- particularly for the case of the moving interface in 
fluid-fluid displacement experiments — still do not fit 
the framework provided by the universality classes we 
discuss. For this reason, we study models in each of the 
universality classes above the depinning transition — i.e., 
in the moving phase. The results obtained allow us to re- 
interpret the experimental measurements. Furthermore, 
we find that for one of the universality classes, some of the 
critical exponents change values at the depinning transi- 
tion. 

The paper is organized as follows: In Sec. II we intro- 
duce the problem, by describing the set of relevant scaling 
properties and associated exponents, for driven interfaces 
moving in disordered media. In Sec. Ill we discuss the 
experimental and numerical results motivated our study. 
In Sec. IV we consider several numerical models, and 
identify two distinct universality classes. In Sec. V we 
try to link the results obtained for the two universality 
classes with the experimental results. To achieve this 
goal, we study models, representative of each of the uni- 
versality classes, above the depinning transition. Finally, 
in Sec. VI, we summarize our results and discuss some 
of the new questions arising from this work. 



II. SCALING PROPERTIES AND CRITICAL 
EXPONENTS 

An interface moving in a disordered medium becomes 
rough due to the action of the disorder. This roughening 
process can be quantified by studying the global interface 
width 



W(L,t) = ( /i 2 (x,t) -/i(x,i) 



1/2 



(2.1) 



where L is the system size, the overbar denotes a spatial 
average, and the angular brackets denote an average over 
realizations of the disorder. 
The study of discrete models 

and continuum 

growth equations [^2|,^3) leads to the observation that 
during the initial period of the growth, i.e. for t<( x (L), 
the global width W grows with time as 



W(t) ~ ^ 



(t<tx), 



(2.2) 



where j3 is the growth exponent. For times much larger 
than t x the width saturates to a constant value. It was 
observed that the saturation width of the interface W sa t 
scales with L as 



(t»t x ), 



(2.3) 



where a is the roughness exponent. 

Th e dep end ence of t x on L allows the combination of 
(|2.2D and (2.3) into a single scaling law pl[ 



W(L,t)~L a f(t/t x ) 



where 



(2.4a) 



(2.4b) 



Here z = a//3 is the dynamical exponent, and f(u) is a 
universal scaling function that grows as u@ when «« 1, 
and approaches a constant when u 1. 

An alternative way of determining the scaling expo- 
nents is to study the l ocal width w in a box of length 
£ < L. The scaling law (^.4j), and the fact that the inter- 
face is self-affine (and with a < 1), allow us to conclude 



w(£,t)~£ a f(£/£ y 



where 



or 



,1/2 



(*«t X ), 



L (t»t x ). 



(2.5a) 



(2.5b) 



(2.5c) 



Here /(it) is a universal scaling function that decreases 
as u~ a when u 3> 1 and approaches a constant when 
u< 1. 

Not all of the exponents introduced so far are indepen- 
dent; e.g., we can calculate the velocity exponent from 
the other scaling exponents, as we show by the following 
argument. Near the depinning transition, except for a 
few growing regions, most of the interface is pinned by 
some pinning path. The growth occurs by the propaga- 
tion of these growing regions through the system. The 
characteristic time required for this propagation must 
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be of the order of t x , since this is the time for correla- 
tions to propagate across the entire system, as defined by 
( |2.4| b). During this process the interface advances from 
one blocking path to the next, the distance advanced is 
typically of the order of W" S a t- Sin ce, clos e to the transi- 
tion, £ ~ L, we can use ( |l.2| ), fl2.3| ), and (2.4) to obtain 



w s&t i t x ~ r / e ~ r 



/ 



— i/(a— 2) 



(2.6) 



Upon comparison with (1.1), we conclude that 

9 = v(z-a). (2.7) 

This exponent "scaling law" can also be derived in 
a different way |4^,|5(|, by considering the interplay be- 
tween the range of action of the quenched disorder and 
the length scales for which the disorder can be consid- 
ered to be time-dependent because of the motion of the 
interface. The quenched disorder is a function of x and h, 
but we can consider the change of variables to a reference 
frame moving with the velocity of the interface 



r)(x, h) — ► rj{x, vt + h). 



(2. 



Now, we can consider the conditions for which vt » h. 
This will happen for the length scale characterizing the 
action of the quenched disorder, £. Thus we have 



v~h/t~ f7f ~ c~ z ~ f v ^ z - a \ 

leading to K%. 



(2.9) 



III. DRIVEN INTERFACES IN DISORDERED 
MEDIA 

A. Experiments 

The most commonly measured exponent in experi- 
ments and simulations is the roughness exponent a. In 
the last few years, several experiments in which quenched 
disorder plays a dominant role have been performed. 

(i) Stokes et al. performed fluid- fluid displacement ex- 
periments, and observed that for some conditions a self- 
affine interface would form, while for other, the growth 
would generate a percolation- like cluster Rj . 

(ii) Rubio et al. measured a for fluid- fluid displace- 
ment experiments [Bj, consistently obtaining 0.73 regard- 
less of the velocity of the interface v. Moreover, they 
observed that the total width of the interface decreased 
with v f§. 

(Hi) Horvath et al. performed similar fluid- fluid dis- 
placement experiment, obtaining a s» 0.81 

(iv) Buldyrev et al. studied the imbibition of coffee in 
pap er towels observing a ~ 0.63 for the pinned interface 
[ p~2| , p~3 Jl6 nj] . Similar experiments were performed by 
Buldyrev et al. for 2 + 1 dimensions leading to q « 0.52 



(v) He et al. conducted fluid- fluid displacements exper- 
iments for different velocities of the interface (ji]] . They 
reported values of a varying from 0.6 for high velocities, 
to 0.9 for low velocities. Furthermore, they also observed 
that the total width of the interface decrease with veloc- 
ity@. 

(vi) Zhang et al. performed flameless paper burning 
experiments JTl[ . They reported a sa 0.71 for the mov- 
ing interface. 

Thus the experiments reveal values of the roughness 
exponent as low as 0.6 and as high as 0.9, leading many 
to question whether one can divide all systems into a 
small number of universality classes (Fig. |2|). 



B. Discrete models 

Several discrete models with quenched disorder have 
been proposed to explain the experimental results. 

(i) Early on, Cieplak and Robbins noticed the impor- 
tance of quenched disorder in many surface phenomena 
p7[ . They studied two models, motivated by fluid- fluid 
displacement experiments. The first, called the fluid in- 
vasion model, is a quite realistic model and calculations 
reveal a w 0.8 in 1 + 1 dimensions |p7| , ^8| . The sec- 
ond, called the random field Ising model (RFIM) fl48[ , 
does not have a self-affine phase in 1 + 1 dimensions, so 
no roughness exponent can be calculated. On the other 
hand, for 2 + 1 dimensions, a self-affine regime exists, and 
a roughness exponent of 0.66 was calculated [|7|2§. 

(ii) Kessler et al. integrated the Edwards- Wilkinson 
(EW) equation with quenched disorder, with the purpose 
of explaining the results for the fluid-fluid displacement 
experiments p9| . They reported values of a that vary 
with the driving force, i.e. with the velocity of the in- 
terface. For low velocities they found a ~ 1, and for 
large velocities, a w 0.5 |p9fl . Furthermore, they found a 
decrease of the total width with the velocity of the inter- 
face. Similar results were obtained earlier by Koplik and 
Levine @. 

(Hi) Buldyrev et al. introduced a discrete model to 
explain their imbibition experiments jl2|,[l3],[l^] . They 
found that in 1 + 1 dimensions the scaling behavior of 
the pinned interface can be obtained exactly by mapping 
the interface onto DP. In higher dimensions they found 
that the interface can be mapped onto a directed surface 
(DS) @. In 1 + 1 dimensions, DP and DS are equiv- 
alent, so this model is referred to as directed percolation 
depinning (DPD) model 

(iv) A similar model was introduced independently by 
Tang and Leschhorn (33|. For the DPD class of models, 
it was observed that a ~ 0.63 for the pinned interface 
01111. 

(v) Parisi also introduced a model for interface rough- 
ening; however, no measurement was made of the rough- 
ness exponent |3l[]. Several authors studied Parisi's 
model and found it to possess problems (see Appendix 
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B). For this reason, Amaral introduced an alternate ver- 
sion of the model, which led to exponents identical to the 
ones obtained for the DPD class of models Q . 

(vi) Dong et al. also integrated numerically the QEW 
equation, and reported a rj 0.97 from the study of the 
height- height correlation function |35| . 

(vii) Leschhorn introduced a model, corresponding es- 
sentially to a discretization of the QEW equation, and 
found a w 1.25 from measurements of the global width. 
The disagreement with Ref. |3j| was later explained by 
Leschhorn and Tang, by showing that the height-height 
correlation function is limited to a scaling with a rough- 
ness exponent smaller or equal to 1 Q . 

(viii) Roux and Hansen studied what amounts to a 
self-organized version of Leschhorn's model finding that 
a w 0.86 for the local width and a w 1.22 for the global 
width ]|{|. Using a similar model, Galluccio and Zhang 
found identical results 40 1. 

(x) Csahok et al. studied the effect of multiplicative 
quenched disorder in the scaling properties of the inter- 
face They found (3 « 0.65 and a « 0.8, in close 
agreement with the experimental results of Horvath et 
al. M. 



C. Continuum models 

On the continuum front, some attempts were made to 
explain the previous results. The first theoretical stud- 
ies to introduce quenched disorder were performed inde- 
pendently by Nattermann et al. }49[] and Narayan and 
Fisher f5C | . Analyzing the Kardar-Parisi-Zhang (KPZ) 



equation [ 23| with quenched disorder, they argued that 
the dependence of the noise on /i(x, t) introduces an in- 
finite hierarchy of nonlinearities. They also noticed that 
in the annealed disorder case the KPZ nonlinearity has a 
kinematic origin. Thus A ~ v, and A should vanish at the 
depinning transition. For these reasons, they argue that 
at the depinning transition, it is sufficient to consider the 
equation |48 5M 



dh 

at 



= F + ?7(x, h) 



\7 2 h, 



(3.1) 



where r/(x, h) represents the quenched disorder. This 
equation was studied by means of the functional renor- 
malization group |l8]-^0| , yielding to first order 



In view of the disagreement between the theoretical 
predictions of the EW equation with quenched disorder 
and the results presented in the previous subsections, 
Csahok et al. numerically integrated the KPZ equation 
with quenched disorder 



— = F + 77(x, h) + V 2 h + X(Vhf 



(3.3) 



The numerical integration of ([3^) for the moving phase 
pl[ yielded exponents in agreement with the calculations 
for the models in the DPD universality class. 



IV. THE DEPINNING TRANSITION 

A. Looking for nonlinearities 

The fact that for (Q we can calculate the exponents, 
and they do not agree with the numerical results obtained 



for (3.2), suggests that the nonlinear term A(V/i) may 



play a crucial role in determining the scaling properties 
of the interface. Unfo rtunately no analytical results are 
available for Eq. (3^3). In this section, we investigate the 
presence of the nonlinear term in several of the models 
discussed in the previous section jl2| . 

If we consider an equation of the KPZ-type, and con- 
sider its average over the spatial coordinates, we obtain 



dh/dt = F + rj + V 2 /i + A(V/i) 2 . 



(4.1) 



Now, if we impose an average tilt to the interface, m 
Vy, and make the transformation 



h(x, t) — ► y(x, t) = h{x, t) + mi, 



(4.2) 



we obtain 



v{m) = dy/dt = F + T] + V 2 h + X(Vh) 2 + \m l . (4.3) 

So, it follows that the average velocity of the interface 
changes with the tilt as |5^| 



e/3, 



(3.2a) 



v(m) = v + Am 2 



(4.4) 



v = 1/(2- a), (3.2b) 

and 

z = 2-2e/9, (3.2c) 

where e = 4 — d. Hence a = 1 in 1 + 1 dimensions, in 
disagreement with many of the experimental and numer- 
ical results. Reference |50| also argued that all but (|3.2|c) 
are exact for all orders in e. 



Thus by varying the tilt m, we can study the presence 
of nonlinear terms in the growth equation and calculate 
the effective coefficient A of a given model j5^] . 

For discrete models the tilt can be easily implemented 
introducing helicoidal boundary conditions. In 1 + 1 
dimensions, we simply impose the condition, /i(0, t) = 
h(L, t) - Lm and h(L + 1, t) = h(l, t) + Lm. In d+ 1 di- 
mensions, we use the same boundary conditions, so that 
the tilt occurs only in one direction. 
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B. The DPD universality class 

1. The DPD-1 model 

We begin by applying the method described above to 
the model introduced by Tang and Leschhorn [Q. In 
1 + 1 dimensions, the interface becomes pinned by a DP 
cluster [[l3 33 1, and the critical dynamics is controlled 
by a divergent correlation length parallel to the inter- 
face £ ~ f~ v , with v k, 1.73. This model, referred to 
as "DPD-1", excludes overhangs, and gives rise to a self- 
affine interface at the depinning transition, with a rough- 
ness exponent a rs 0.63. 

In the DPD-1 model, in 1 + 1 dimensions, we start from 
a horizontal interface at the bottom edge of a lattice of 
size L. At every site of the lattice we define a random, 
uncorrelated quenched variable, the noise to with magni- 
tude in the range [0, 1]. During the time evolution of the 
interface, we choose one of the L columns at random |53) . 
If the difference in height to the lowest neighbor is larger 
than +1, this lowest neighboring column grows by one 
unit. Otherwise, the chosen column grows one unit pro- 
vided the noise on the site above the interface is smaller 
than the driving force F. The unit time is defined to be 
L growth attempts. 

We measure the velocity of the interface for different 
values of the reduced force / and different values of the 
tilt to. The results for 1 + 1 dimensions are shown in Fig. 
J(a). For a fixed force /, we find that the interface ve- 
locity depends on the tilt to, indicating the existence of 
nonlinear terms. Near the depinning transition (f — > 0), 
the velocity curves become "steeper" and, from (4.4), we 
infer that A is coupled to the external force and that it 
increases as / — > 0. 

To measure A, we first attempt to fit a parabola to the 
tilt-dependent velocities in the vicinity of zero tilt. The 
calculations indicate that as we approach the depinning 
transition, A diverges as |M2| 



(4.5) 



However, in the vicinity of F c , the velocity curves lose 
their parabolic shape for large tilts (see Figs. |^(a) and 
| i|)^ indicating a change to a different universality class 

We can understand the breakdown of (4.4) for large 



to using scaling arguments. Substituting Eqs. (1.1) and 
(O) into (O), we find E§ 



v(m,f) cx f +af-' t >m 2 



(4.6) 



Equation ([Oj) indicates that the velocity curves corre- 
sponding to two different forces f\ and /2, with f\ > fo, 
will intersect at a tilt m x (see Fig. [|). For tilts greater 
than to x , v(m, fx) < v(m,f2), a clearly unphysical pre- 
diction since /1 > /2, and since the average velocity, for 
the same tilt m, should be larger for the larger force. 
Thus the velocity cannot follow a parabola for arbitrar- 
ily large m, and a crossover to a different behavior than 



that of Eq. ( f4.6| ) must occur for values of the tilt larger 
than to x . 

Letting (fi — ^2) —* 0, we find from (4.6) that the 



crossing point of the two corresponding parabolas scales 
as 



in" 



f 



e+4> 



(4.7) 



Equations ( |4.6|) a nd ( [4.7| ) motivate the scaling form for 
the velocities |U2f 



v(m,f)~fg( m yf+% 



(4.8) 



where g(x) ~ const for x <C 1, and g(x) ~ x e ^ e+ ^ for 
x 3> 1 [pifl . Figure || shows the data collapse we obtain 



using (4.8), with the exponents 

61 = 0.64 + 0.08, ^ = 0.64 ±0.08, 



(4.9) 



for 1 + 1 dimensions. Below, we will show that an entire 
class of models produces similar results, which will allow 
us to group them into a single universality class. 



2. Other models 



The scaling behavior fl4.8| ) is not limited to the DPD-1 
model in 1 + 1 dimensions. For 2 + 1 dimensions and for 
the models introduced in Refs. 12 3^j3^], we find a very 
similar behavior (Table |) . 

The "DPD-2" model was introduced by Buldyrev et 
al. p^-p] (see also Amaral et al. |[7]]). Let us define 
the model for 1 + 1 dimensions. In a square lattice of 
edge L, assign an uncorrelated random number, the dis- 
order rji, with magnitude uniformly distributed in the 
interval [0,1], to each cell i. We compare the random 
pinning forces to in the lattice with the driving force F, 
where < F < 1. If the pinning force at a certain 
cell to is larger than the driving force, the cell is labeled 
"blocked" , otherwise it is labeled "unblocked" . Thus a 
cell is blocked with a probability p = 1 — F. 

Since the model was developed to study imbibition, we 
will refer to the growing, invading, region as "wet" , and 
to the invaded region as "dry". At time t = 0, we wet all 
cells in the bottom row of the lattice. Then, we select a 
column at random [ p3| , and wet all dry unblocked cells in 
that column that are nearest-neighbors to a wet cell. To 
obtain a single-valued interface, we impose the auxiliary 
rule that all dry blocked cells below a wet cell become wet 
as well fss]] . We refer to this rule as erosion of overhangs. 
The time unit is defined as L growth attempts. 

Figure [| shows the dependence of the velocity on the 
tilt for the "DPD- 2" m odel and the respective data col- 
lapse according to ( |4.8| ), for 1 + 1 dimensions. The values 
of the exponents in the data collapse are |5q] 



= 0.59 + 0.12, 



= 0.55 + 0.12. 



(4.10) 
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Finally, we simulate the model proposed by Amaral 
Hi, which we will refer to as "DPD-3". We simulate it 
in 1 + 1 dimension s, an d rescale the tilt dependent_yeloc- 
ities according to (4. 



6 = 0.70 ±0.12 



using the exponents (Fig. |^) 

= 0.65 ±0.12. (4.11) 



forward growth. Furthermore, Makse argued that the 
application of constraints to the local slopes of the inter- 
faces, as is the case for the DPD class of models, appears 
to be an essential ingredient for A to diverge |5jJ. 



C. The QEW universality class 



3. Discussion 



1. Discretizations of the QEW equation 



Three questions are raised by these results: 

(i) Is the new exponent <f> an independent exponent or 
can it be expressed in terms of the other known expo- 
nents? 

(ii) What are the scaling properties of the tilted inter- 
face? 

(in) What are the reasons for the divergence of A at 
the depinning transition? 

Questions (i) and (ii) have recently been answered by 
Tang, Kardar and Dhar (TKD) Q. They noticed that 
for the DPD class of models there is a characteristic slope 



v{l-a) 



(4.12) 



where £j_ and £|| are, respectively, the perpendicular and 
parallel correlation lengths of a DP cluster (see Appendix 
A)- 

For values of the average tilt smaller than m c , the inter- 
face still belongs to the DPD universality class. However, 
for larger values of the tilt a qualitative change should oc- 
cur. But we showed above that such a change occurs for 
m > nix, so we conclude that m c = m%, from which 
follows Ii 



6 + <p = 2z/(l - a). (4.13) 

Combining fl4.13| ) with (|TTT|) , we find the scaling law 

<f> = v(2 -a- z). (4.14) 

TKD also showed that for m > m x a new universal- 
ity class is obtained Q| , and they were able to relate its 
scaling exponents to the exponents of the KPZ equation 
with annealed disorder |43l 



Otd 



KPZ /KPZ 
y d-l I z d-l 



Zd 



(4.15) 



To address question (Hi), let us consider the effect of 
the nonlinear term. Its role is to make the interface grow 
wherever it has a nonzero slope, so it corresponds to 
what we can refer to as lateral growth. When an in- 
terface is moving, we can consider two types of growth, 
forward growth and lateral growth. The meaning of 
A — > oo is then the irrelevance of forward growth com- 
pared with lateral growth when we approach the depin- 
ning transition. These ideas were studied quantitatively 
by Makse fl57| , who found that, near the depinning transi- 
tion, the probability of occurrence of lateral growth van- 
ishes much slower than the probability of occurrence of 



We mentioned in Sec. Ill that several authors numer- 
ically integrated Eq. ( |3.lD , and others developed mod- 
els corresponding to the discretization of this equation. 
Thus, we expect results in agreement with the theoreti- 
cal calculations. Furthermore, we do not expect, for these 
models, to detect any dependence of the velocity on the 
tilt of the interface. In fact, we have simulated the model 
introduced by Leschhorn f34j, which we will refer to as 
"QEW-1", and find that for any reduced force, A = 
(Fig. @). 

A Hamiltonian model, which we will refer to as "QEW- 
2", was introduced by Makse |jl6). In 1 + 1 dimensions, 
the Hamiltonian is defined as 



n = ]T [{h i+ i - h t f - Fhi + ((i, h)] 



(4.16) 



Here, the first term represents the elastic energy that 
tends to smooth the interface, and hi) is an uncor- 
rected random number which mimics a random poten- 
tial due to the disorder of the medium. In the simula- 
tion, a column i is chosen and its height is updated to 
(hi + 1) if the change in ( 4.16 ) is negative. Thus, only 
motions that decrease the total energy of the system are 
accepted. Backward motions are neglected since these 
are rare events. Equation (3.1) can be easily derived 
from this Hamiltonian simply by considering that 



dh 
dt 



5H 
lh 



a j- si, ' r ? si, ' 



5h' 



(4.17) 



We study the dependence of the velocity with the tilt 
for this model and find that A = for any force, just as 
for the QEW-1 model. 



2. The RFIM 

An interesting result, presented in Sec. Ill, is the rea- 
sonable agreement between the exponents measured for 
the FIM and the RFIM with the predictions of 

( |3 . l[ ) - Since these models are not simple discretizations 
of Eq. (3.1), it is interesting to investigate the presence 
of any nonlinear term. 

For this reason, we study the RFIM in both 1 + 1 and 
2 + 1 dimensions. This model allows for overhangs, and 
for certain values of its parameters can be mapped to the 
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isotropic percolation problem pq| . In the RFIM, spins 
on a square lattice interact through the Hamiltonian 



n = -J2 s i s j-J2i F + ^ i ' h ^ 



(4.18) 



{id) 



where Si = ±1, F now denotes the external magnetic 
field, and £ is the time-independent local random field 
(i.e., quenched noise) whose values are uniformly dis- 
tributed in the interval [—A, A]. The strength of the 
quenched disorder is characterized by the parameter A. 
At time zero, all spins are "down" — except those in the 
first row, which are initially "up". The interface con- 
sists of all down spins that are nearest-neighbors to an 
up spin. During the time evolution of the system, we flip 
any down spin that belongs to the interface and is "un- 
stable," i.e., whenever the flip will lower the total energy 
of the system. The control parameter of the depinning 
transition is the external magnetic field F; the unit time 
corresponds to the parallel flipping of all unstable spins 

For 1 + 1 dimensions, there are two morphologically- 
different regimes, depending on the strength A of the 
disorder (i.e., of the random fields). For A > 1.0, the in- 
terface is self-similar (SS), while for A < 1.0 it is faceted 
or flat (FA). For 2 + 1 dimensions, there is again a FA 
regime (A < 2.4), a SS regime (A > 3.4), and also a self- 
afline (SA) regime in between (2.4 < A < 3.4) The 
SA regime, which exists only for the case of 2 + 1 dimen- 
sions, is t he o nly regime of the RFIM for which either 
Eqs. (3.1) or (3.3) could apply. In the SS regime, over- 
hangs in the interface are larger enough for the approxi- 
mation of small slopes not to be valid. On the other hand, 
in the FA regime, lattice effects dominate the grow th [ |60[ . 

For the SA and SS regimes we find that (fD|) is still 
valid; however, we find a negative <f> 

A~/W-f0. (4.19) 

This behavior can be understood, for the SS regime, by 
considering that near the depinning transition, the mor- 
phology of the interface corresponds to the hull of an 
isotropic percolation cluster, which has no well-defined 
orientation p7[ . Thus a change in the boundary condi- 
tions will not affect the growth process, and we cannot 
expect any divergence of a possible nonlinear term when 
the magnetic field approaches its critical value. On the 
other hand, for large fields the effect of the quenched 
disorder diminishes, and we can observe an average in- 
terface orientation. For such values of field, we expect 
the presence of nonlinear terms, of kinematic origin, to 
be relevant. Although for the SA regime the behavior of 
A is similar (Fig. ||), we do not have a simple argument 
that would explain the observed behavior. 

These results lead us to conclude that in the SA regime 
the RFIM belongs to the universality class of Eq. (3.1). 
This conclusion is further supported by the agreement 
between the numerically determined exponents, a ~ 0.67 
and 9 rj 0.60 for 2 + 1 dimensions, and the theoretical 



D. Discussion 



The results of Table | show, for 1 + 1 dimensions, a 
separation into two groups in the values of the critical 
exponents for the six models studied (5(|. This sepa- 
ration reflects the existence of two distinct universality 
classes, described by the two continuum growth equa- 
tions, fO) and (|OD. For the QEW-1, QEW-2 models, 



and for the RFIM in the SA regime, we find that ei- 
ther A = or A — ► at the depinning transition (Fig. 
|9|). Thus, the scaling behavior of these models should 
be correctly described by (Pi) . For the DPD-1, DPD-2, 
and DPD-3 models we observe a divergent A, indicating 
that nonlinearities are relevant near the depinning tran- 
sition. Thus to properly describe the sc alin g prope rtie s 
of these models it is necessary to study (|3~3|), since (3.1) 
does not include the nonlinear term A(V/i) 2 . Further 
evidence of the existence of the two universality classes 
is given by the values of the roughness exponents. The 
models for which A diverges at the depinning transition, 
have a ks 0.63, in agreement with the mapping to DP 
p9] |. O n the other hand, models in the universality class 
of Eq. (3.1), gave roughness exponents typically larger, 
and in better agreement with the theoretical prediction, 
awl. 

It is worthwhile to discuss the implications of the pre- 
vious results on other models and theoretical approaches 
studied in the literature. Recently, considerable atten- 
tion has been focused on a self-organized version of 
the DPD model, the so-called self-organized depinning 
(SOD) model @,0,|l]-|6|] . In this model, the constant 
driving force of the DPD model is replaced by an algo- 
rithm that selects the site on the interface for which the 
disorder is weakest — i.e, it selects the point in the inter- 
face offering the least resistance. Thus, the local update 
rule of the DPD model is replaced by a global search pro- 
cess, similar to the one used in invasion percolation. This 
update rule assures that the model will evolve to criti- 
cality and that it will remain in the critical state. This 
implies that the SOD model is always at the depinning 
transition. Thus the properties of the SOD model are de- 
scribed by the DPD model at criticality, which implies, 
since the SOD model is always at criticality, that the 
SOD model cannot be described by a continuum equa- 
tion of the KPZ-type, since this would require an infinite 
coefficient of the nonlinear term. However, this argument 
does not rules out that the SOD model can be described 
by a continuum equation that upon a mapping to the 
quenched KPZ equation will generate an effective coeffi- 
cient A that would diverge [|43|. 



V. THE MOVING PHASE 



predictions of (3.1 
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A. Why look at the moving phase? 



<PA = <Pm + v(a M - a A ). 



(5.4) 



In the previous section we identified the universal- 
ity classes for interface roughening in the presence of 
quenched disorder. However, the experimental situation 
still leaves us with several unsolved puzzles, such as the 
change in the roughness exponent or the dependence of 
the global width on the driving force for the fluid-fluid 
displacement experiments. 

In trying to shed some light over the remaining prob- 
lems, we observe that most of the numerical and the an- 
alytical studies focus on the "pinned phase" [F < F c ) 
while nearly all experimental results are for the "mov- 
ing phase" (F > F c ). A logical next step is therefore to 
study the moving phase for models representative of each 
of the two universality classes. 



B. The DPD universality class 

1. The steady state 

Let us start by considering the DPD universality class. 
In Fig. |l^(a), we show the local width for the pinned and 
the moving phases. We find otp sa 0.63 for the pinned 
phase, and ccm w 0.75 for the moving phase (Table |J). 
So far, no explanation for this change in roughness expo- 
nent at the depinning transition is available. 

The study of the scaling of the local width for the mov- 
ing phase reveals a novel dependence of the prefactor 
of the width on the external force. We also find that 
the exponents characterizing this anomalous dependence 
change values for the two different regimes. In the first 



w 



f 



-<PM 



{t « 0: 



(5.1) 



where £ is the correlation length, and tpM is a new expo- 
nent that characterizes the dependence of the prefactor 
of the width on the driving force. In the second regime, 
the effect of the quenched disorder becomes irrelevant 
compared to the annealed disorder |53|| , and we obtain 



w 



f 



-<PA 



(5.2) 



where ifiA is a new exponent, and a a is the roughness ex- 
ponent corresponding to annealed disorder. Depending 
on the absence or presence of nonlinear terms, we recover 
the results of either the EW or the KPZ equations with 
annealed disorder. This crossover was observed both for 
experiments and simulations of discrete models |5j . 

Due to the dependence of the prefactor of the width 
on the driving force, and since £ ~ , we propose the 
scaling Ansatz 



(5.3) 



Upon comparison with (5.1) and (5.2), we find that the 
scaling function g(u) satisfies g(u 3> 1) ~ const, and 
g(u -C 1) ~ u aM ~ aA . Using (5.3) we also obtain 



In Fig. [n](b) we show the data collapse of the widths for 
different forces and length scales. The deviations from 
scaling for large values of t/^ are due to finite-size effects 
(see Fig. |l0|(d) for discussion). 

While previous studies considered clm as an effective 
exponent, and concluded that no self-affine scaling exists 
in the moving phase for i <c £ (3^], we find a consistent 
value for clm regardless of how closely we approach the 
depinning transition. The self-affinity of the interface is 
supported by the good data collapse obtained with flS.j), 
as shown in Fig. nOl 

To derive a second relation for the new exponent tfM, 
let us consider the approach to the depinning transition 
for a system of size L. For a finite system, the transition 
does not occur for / = 0, as it would for an infinite sys- 
tem, but rather for an effective critical force F C {L) such 
that i ~ \F C {L) -F c \- V ~ L, implying / ~ L^". Using 
this result and (|5.l|), we obtain 



(5.5) 



For the pinned phase we have 

w~l ap (£<L~£). (5.6) 



At the depinning transition, ( |5.5| ) and (5.6) should scale 
in the same way. If we replace I by L we obtain Wm 



<Pm = v{a P - at M )- 



(5.7) 



Rep lacing the known values of ap, cum, and a a into 
(|j) and (|3), we find 



-0.20, if A w 0.23, 



(5.8) 



to be compared with the measured values ip&i w —0.12 
and (fiA ~ 0.34 (Table |l|). Although the agreement with 
the measured values is not p erfect, the error bars do not 
rule out the validity of (|5.7|). 



2. The time evolution 

For the DPD universality class, we find that (3 changes 
value at the depinning transition (see Table 0). We find 
[3p w 0.63 for the pinned phase and Pm ~ 0.75 for the 
moving phase. We obtain for both sides of the transi- 
tion a k, (3 implying that z remains unchanged at the 
transition |6q| . The exponent z characterizes the time 
scale t x for th e pr opagation of correlations in the inter- 
face (cf. Eq. (2.4)). This time scale is not expected to 
depend on the external force p4| . These results a re in 
good agreement with a numerical integration of Eq. ( |3.3| ) 

13- 

For the time evolution of the width, we again find an 
anomalous dependence of the width on the reduced for 
the two different scaling regimes. In regime I, 
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w^t >3M r KM (<<*x), 



(5.9) 



where t x ~ £ 2 is the characteristic time for the prop- 
agation of correlations over a distance £, and km is a 
new exponent that characterizes the dependence of the 
prefactor of the width on the driving force. In regime 
II, the effect of the pinned disorder becomes irrelevant 
compared to the annealed disorder |53[] , and we obtain 



W~tP A f-" A (t»t x ), 



(5.10) 



where ka is a new exponent, and Pa is the growth expo- 
nent corresponding to annealed disorder. 

Since t y ~ ~ we propose the scaling Ansatz 



w(t,f)~tf>* r KA 9(t/tx). 



(5.11) 



Upon comparison with ( [5.9] ) and ( 5.10| ), we find that the 
scaling function g(u) satisfi es q( u 3> 1) ~ const, and 
g(u Cl)~ u 0M ~P A . From (|5.1l[) we also obtain 



ka = «m + 2^(/3 M - Pa)- 



(5.12) 



In Fig. [ll] we display the rescaling of our simulation re- 
sults for the dependence of the g lobal width on time and 
driving force according to ( 5.11 ). 

To determine a second relation for the new exponent 
Km, let us again consider the approach to the depinning 
transition for a system of size L. As discussed above, the 
transition occurs for an effective critical force such that 
£ ~ L, im plying / ~ L^ 1 ^ ~ t^^ zh ' . Using this result 
and (5.9), we obtain 



For the pinned phase, we have 

W~t 0p (t<t x ). 



(5.13) 



(5.14) 



At the transition, fl5.13 ) and ( 5.14 ) must scale in a similar 
way. So, if we replace t by t x , we obtain 



KM = ZV>(Pp - P M )- 



(5.15) 



Repla cing t he kn own values of v, z, Pp, Pm, and Pa into 
(PJ) and (|5~T5|), we find 



km 



-0.12, K A w 0.64, 



(5.16) 



in good agreement with the simulation results, km 
-0.11 and K A w 0.65 (Table [fi|). 



local width on the system size L, that can be described 
by the scaling function [p6pql 



(5.17) 



w(t,L) ~ L aG $(£/L), 



where $(u) is a scaling function that for «<1 scales as 
u ap , and ac is the global roughness exponent, which is 
defined by the scaling of the gl obal width with the sys- 
tem size: W ~ L a ° . Equation ( 5.17 ) is valid only when 
a G > 1, and it implies that w(£,L) ~ f p L aG -* p for 
l«L. Hence, the local slopes are unbounded, i.e., they 
diverge with the system size. This anomalous charac- 
teristic of the QEW models is displayed in Fig. p"2|(a), 
where we plot interfaces generated for the QEW-2 model 
for different values of the system size. We calculate the 
global width as a function of L and find at the depinning 
transition that ag ri 1.23. 

We also determine the scaling of the local width w in a 
window of size £, for different values of the reduced force. 
In Fig. |l3|(a), we show the local width for the pinned and 
the moving phases. Analysis of consecutive slopes yield 
roughness exponents ap rj 0.92 for the pinned phase, 
and o>m ~ 0.92 for the moving phase (Table ||). These 
results are in agreement with values commonly found for 
the QEW class of models when the roughness exponent 
is calculated from the local width — i.e., ap rs olm ~ 1- 
The slightly smaller value of our estimates is is likely due 
to finite size effects. 

As for the DPD universality class, the study of the local 
width in the moving phase reveals an anomalou s de pen- 
dence on the reduced force, which leads to Eqs. ( |5.3|) and 
( |5.4| ) being also valid for the QEW universality class. In 
Fig. [l^(b), we s how the data collapse of the local widths 
obtained using (5.3). 

To determine a second relation for the new exponent 
ifM, we again consider the approach to the depin ning 
transition of a system of size L. As before, we obtain ( |5.5| ) 
for the moving phase. However, for the pinned phase, the 
scaling behavior (5.17), leads to 



gap jjy.Q—ap 



(£<^L). 



(5.18) 



At the transition, ( |5.5|) and ( [5.18 ) must be identical. 
Hence, we find 



and 



ap = aM, 



ip M = v(a G - a P ). 



(5.19a) 



(5.19b) 



C. The QEW universality class 

1. The steady state 

For the study of the QEW universality class, we focus 
on the QEW-2 model. We observe a dependence of the 



Rep lacin g the mea sured values of is, a G , ap, and a a 
into $% and fl5.19|), we find 



ip M w 0.42, ip a » 1.04, 



(5.20) 



in good agreement with the values obtained in our simu- 
lations, ip M w 0.44 and tp A « 0.99 (Table ||. 
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2. The time evolution 

We study the scaling of W with t, and find (3p 0.85 
and Pm ~ 0.86. These results imply that z remains un- 
changed at the transition. In Fig. [14], we show the scal- 
ing of the global width as a function of time for different 
value s of the reduced force, and its data collapse using 



For the QEW universality class, we find by repl acing 
the measured values of v, z, [3p, and [3a into (5.12) and 



(|5.15[) that 



km « 0.02, k a «1.19, 



(5.21) 



in good agreement with the values obtained in our simu- 
lations, km ~ and ka ~ 1-15 (Table ||). 



D. Discussion 

The results obtained in the previous subsections al- 
low us to discuss the experimental and numerical results 
reported in the literature. A problem with the inter- 
pretation of experimental and numerical results for the 
roughness exponent, has been the wide range of values for 
a: 0.5 — 1.25, measured in the moving phase. We note 
that, in this regime, the crossover to the annealed dis- 
order regime leads to "effective" exponents that change 
with the velocity (or the driving for ce). For this reason, 
we suggest that the scaling function ( |5.3|) might be useful 
in the determination of the exponents from the study of 
the local width w. As shown in Table [n], the exponent 
ifM has different signs for the two universality classes, 
leading to sharply distinct scaling behaviors for the pref- 
actor of the width with the driving force. Since in many 
experiments it is possible to monitor the velocity of the 
interface, and therefore the driving force, the study of 
this prefactor may also lead to an easier identification of 
the universality class to which the experiments belong. 

In light of this discussion, the interpretation of the re- 
sults ofRef. 1 29 1 is clear. The numerical integration of the 
EW equation with quenched disorder performed in Ref. 
|p9| must belong to the QEW universality class, and the 
reported exponents are "effective" exponents whose val- 
ues were affected by the crossover to the annealed regime. 

The determination of the universality class of the fluid 
invasion model (FIM) of Ref. |27],|28| is not trivial. How- 
ever, we note that v is identical for the FIM and for the 
Hamiltonian model. On the other hand, the value of 
a measured in Ref. 0,^8| is somewhat smaller than the 
value obtained for the QEW model using the local width. 

Regarding the fluid-fluid displacement experiments of 
Refs. f7[-^0|, we find that they share several scaling prop- 
erties with the Hamiltonian model. We first notice that 
the range of roughness exponents measured in the exper- 
iments of Refs. [f7| [l0|] is consistent with the values that 
could be obtained with the Hamiltonian model. More- 
over, Rubio et al. found that the amplitude of 



the width decreases with the increase of the velocity. 
This property is entirely consistent with our findings for 
t he QEW universality class: if we replace / ~ v 1 / 8 in 
( |5.l|) , then we obtain that the width decreases with v as 

W ~ V ~Vm/S _ 

The imbibition experiments of Refs . , |l3| , [il| and the 
paper burning experiments of Ref. [|l lfare believed to 
belong to the DPD universality class. This conclusion 
is supported by the values of the exponents measured in 
both experiments. 

The reason for the different universality classes for the 
two groups of experiments, fluid-fluid displacement on 
one hand and imbibition and paper burning on the other, 
is not totally clear. However, a possible explanation 
might be the importance of lateral growth in the second 
case compared with the first. When we analyze the imbi- 
bition and paper burning experiments, it becomes clear 
that the interface can be described by a single- valued in- 
terface, because only the highest position of the interface 
in each column is relevant for the dynamics. This ef- 
fectively corresponds to eroding any overhangs and, as 
shown by Makse J57J , leads to a diverging A at the tran- 
sition. On the other hand, in the fluid-fluid displacement 
experiments, we must at all times consider the full inter- 
face, because of the trapped fluid that might exist below 
the interface. Since it will require some effort to displace 
this fluid or to separate it from the rest of the displaced 
fluid, we cannot physically justify a rule such as the ero- 
sion of overhangs. This implies that A vanishes at the 
depinning transition, and the experiments must belong 
to the QEW universality class. 



VI. SUMMARY 

To summarize, we perform a systematic study of sev- 
eral models proposed to understand the motion of driven 
interfaces in disordered media. We are able to identify 
two distinct universality classes. For one of these univer- 
sality classes, QEW, we observe that it can be described 
by a Langevin equation similar to the EW equation but 
with quenched disorder. For the QEW universality class 
nearly all exponents for the pinned phase can be obtained 
by renormalization group calculations. 

For the other universality class, DPD, we find that 
it can be described by a Langevin equation similar to 
the KPZ equation, but with quenched disorder. Further- 
more, we find that the coefficient of the nonlinear term 
A diverges at the depinning transition. Because of this 
divergence, no perturbative analytical calculations have 
been possible so far. However, a mapping of the static 
properties of the interface, at the depinning transition, 
to DP yields all static exponents of the problem. A map- 
ping of the dynamics to isotropic percolation yields the 
dynamical exponent. 

The identification of these two universality classes 
leads to some understanding of many of the results ob- 
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tained previously for models and experiments. However, 
many other results could not be understood in terms of 
the exponents obtained for the two universality classes 
at the depinning transition. For this reason, we also in- 
vestigated the scaling properties of the models in the two 
universality classes above the depinning transition. 

We find that for the DPD universality class, a and j3 
change their values at the depinning transition. Further- 
more, we find a dependence of the prefactor of the width 
on the driving force. The new exponent y>Mi character- 
izing the scaling of this prefactor, can be used to relate 
the values of the roughness exponents on both sides of 
the transition. 

For the QEW universality class, we find that a and (3 
remain unchanged at the depinning transition. As for the 
DPD universality class, we also find a dependence of the 
prefactor of the width on the driving force. The exponent 
ipM is, in this case, found to relate the different values 
of the local, ap, and global, ac, roughness exponents at 
the depinning transition. 

These results provide us with a more consistent under- 
standing of the scaling properties of the two universality 
classes both at and above the depinning transition. This 
knowledge enables us to interpret most of the experimen- 
tal results obtained previously. However, several ques- 
tions, that demand an answer, are still unsolved. One of 
them is the reason for the change in values of the rough- 
ness exponent at the depinning transition for the DPD 
universality class. Another is the understanding of the 
different values of ac and ap for the QEW universality 
class. 
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APPENDIX A: PREDICTING THE EXPONENTS 
OF THE DPD MODEL 

1. Directed percolation 

Near p c , the size of DP clusters is characterized by a 
longitudinal (parallel) correlation length £|| and a trans- 
verse (perpendicular) correlation length £j_ that diverge 
as 0,0 



The parallel and perpendicular correlation length expo- 
nents for DP clusters have been calculated |7^] , with the 
results 



1.733 ±0.001, t/i = 1.097 ±0.001 



(d=l). 
(A2) 



2. Static properties 



The mapping of the pinned interface to DP enables us 
to estimate the static exponents of this problem from the 
characteristic exponents of DP clusters. The character- 
istic length £ of the pinned regions must be of the order 
of £|| , so we can identify the exponent v to be 

v = v v (A3) 

The global width W sa t of the pinned interface should 
scale as since its advance is blocked by a DP path. 
On the other hand, £h must be larger than the system 
size L for the interface to become pinned, from which 

follows ITlilSiil 



W sl 



"a /"II 

II 



L VJ J"* > L). (A4) 



Comparing with (2.3), we conclude that the roughness 
exponent is given in terms of the correlation exponents 
for DP, 

a = v±/vu. (A5) 
Substituting (|a|) into ([A5|), we predict 

a = 0.633 ±0.001 (d=l). (A6) 

3. The minimum path 



For any fractal set with a fractal dimension smaller 
than the dimension of the embedding space, the distance 
between to points of the set I is not given by the Eu- 
clidean distance r between those points but by the mini- 
mum path distance (or chemical distance) l m j n . Numer- 
ical studies show that £ m -, n scales as 



(A7) 



€\\ ~ \pc-p\ 



£-L ~ \Pc ~P\ 



(Al) where d m i n is called the minimum path exponent. 
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4. Dynamical properties 

Recently, Havlin et al. showed that the dynamics of 
the DPD model in d + 1 dimensions could be mapped 
to the minimum path of isotropic percolation in d di- 
mensions [Q. They argued that close to the depinning 
transition only a few cells on the interface are active, i.e., 
unblocked. These active cells perform a global search of 
unblocked cells and propagate correlations in the inter- 
face. Since the cells through which the active cells can 
move are confined to unblocked or eroded blocked cells 
in a region of thickness £j_ -C £11 , we can identify this 
region with an isotropic percolation cluster embedded in 
a space of dimension d. 

Thus, the path through which active cells move is the 
minimum path. So, the time t for the correlations to 
propagate an Euclidean distance r, is proportional to 
Anin, which scales as r dmin for isotropic percolation. From 
this follows that we can connect the dynamical exponent 
z to the static exponent d m i n of isotropic percolation 

z = d min . (A8) 



APPENDIX B: THE DPD-3 MODEL 

In this appendix, we discuss the model introduced by 
Parisi to study interface roughening in the presence of 
quenched disorder ]3l| . We show that the original growth 
rule leads to an unphysical behavior. We then proceed to 
describe the model proposed by Amaral that solves the 
problems with the original model |32| . 

Let us start by describing the model originally intro- 
duced by Parisi. In a square lattice of edge L, with pe- 
riodic boundary conditions in that direction, we assign 
to each cell an uncorrelated random number, uniformly 
distributed in the interval [—1,1]. We start, at t = 0, 
by invading the bottom row in the lattice, and then, we 
calculate the driving force for each column 



F(i, t) = p max[/i(j ± 1, t) + 1 - h(i, t)}, 



(Bl) 



where p is the controlling parameter of the model. If this 
driving force Fi is larger than the noise in the cell above 
the interface, r](i, h(i) + 1), then we update the height of 
column i to 



is difficult to justify. In fact, it does not seem reasonable 
that any growth mechanism would pull a lower column 
to an height larger than that of its neighbors. 

For this rea son, Amaral proposed a new growth rule 
to replace (B2) 



h(i,t + l) = h(i,t) + 1. 



(B3) 



In this way, the invading region only advances one cell at 
a time; thus, avoiding the formation of big jumps between 
neighboring columns |32|. 

The study of this model, which we will refer to as 
"DPD-3" , leads to the result that a critical value of the 
controlling parameter exists, p c = 0.196 ± 0.002. Be- 
low this threshold the interface stops moving after some 
finite time, and above it it moves indefinitely. To de- 
termine some of the scaling exponents, we applied the 
gradient method, introduced by Sapoval et al. Jr3 |, and 



already used in the study of the DPD-2 model [10 17 
and in several other studies j74|-|76| . 

Figure [l|| displays our results for the case of a linear 
gradient, Vp(h) = g. The scaling of the interface changes 
with the concentration rate and is characterized by an 
exponent 7 according to 



w(l,g)~£ a f(e/g-<^ 



(B4) 



where / is an universal scaling function that satisfies 
f(x < 1) ~ const, and f(x > 1) - x~ a . Figure |l|(a) 
shows the width w as a function of the size £ for differ- 
ent values of the concentration gradient g. Figure |l5|(b ) 
shows the data collapse of these results according to 1 
The best collapse was obtained with the exponents 



a = 0.60 ±0.05, 7 = 0.49 ±0.05. 
Since we know that [lira] 



a = u±/un, 7 = f j_/(l + v±), 



(B5) 



(B6) 



we obtain 



z/_l = 1.0 ± 0.1, v\\ = 1.6 ±0.1, (B7) 
in agreement with the known results for DP 173] . 



h(i,t+ 1) = m-Ax[h(i ± l,i) + 1]. 



(B2) 



The time evolution of this model leads to the develop- 
ment of huge jumps in the height of neighboring columns 
that propagate along the system. For small enough val- 
ues of p, only one or two such growing regions will be 
propagating in the system. The reason for suc h "u nnat- 
ural" behavior of the model is the growth rule ([B2|). The 
fact that the advancing column grows to an higher value 
than the highest neighbor is a kind of "bootstrap" that 
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FIG. 1. The depinning transition. In the "pinned phase," 
F < F c , the velocity of the interface is zero. In the "moving 
phase," F > F c , the interface moves with an average velocity 
v = v(F), where v(F) ~ (F — F c ) e for F close to F c , and 
v(f) ~ F for F F c . Thus, the velocity plays the role of the 
order parameter of the transition. 



14 



FIG. 2. Values of the roughness exponent a from most of 
the experiments,simulations and theoretical works reported in 
the literature. The experimental points are, reading from left 
to right, from Rcfs. [8], [9], [10] (a range of values is reported 
of which the extremes are plotted), [12-17], [11]. The simula- 
tion points are, reading from left to right, from Refs. [17] (a 
range of values is reported of which the extremes are plotted) , 
[29], [12-17], [32], [31], [34], [33], [40], [27-28], [45], [46]. Fi- 
nally, the theoretical points, again reading from left to right, 
are from Refs. [48] and [49]. Visually apparent is the wide 
spread of the results. To guide the eye, we indicate the values 
of a predicted by the KPZ equation with annealed disorder 
(bottom line) , the DPD universality class in the pinned phase 
and moving phases (middle lines) , and the QEW universality 
class (top line) . We can see that most results for simulations 
and theoretical calculations are close to either the DPD or the 
QEW line. However, for the experimental results, the agree- 
ment is not so good. In the text we argue that the reasons for 
these disagreements are the crossover present in the moving 
phase for both universality classes, and the fact that for the 
DPD universality class am =fc «p. 

FIG. 3. DPD-1 model, (a) Dependence on the tilt m of 
the average velocity in 1 + 1 dimensions. Data for different 
values of the force / = (F — F c )/F c are indicated by differ- 
ent symbols, ranging from 0.016 (bottom curve) to 0.350 (top 
curve). The system size is L = 512 and each result was av- 
eraged over 30 realizations of the disorder, (b) Data collapse 
of the velocities according to (4.7) using the values for the 
DPD-3 model exponents (Table I). 

FIG. 4. Here we exemplify the "noncrossing" effect on the 
velocity parabolas. We show a perfect parabolic behavior for 
two different forces, /i > /2 (dashed lines) as predicted by 
Eq. (4.3). Also shown is the "curving back" of the velocity 
curve for the smaller force fi (solid line) which is necessary 
in order not to cross the velocity curve for fi . 

FIG. 5. DPD-2 model, (a) Dependence on the tilt m 
of the average velocity in 1 + 1 dimensions. Data for differ- 
ent forces / are indicated by different symbols, ranging from 
0.0149 (bottom curve) to 0.0719 (top curve). The system size 
is 512 and each result was averaged over 30 realizations of the 
disorder, (b) Data collapse of the velocities according to (4.7) 
using the values for the DPD-2 model exponents (Table I). 



FIG. 6. DPD-3 model, (a) Dependence on the tilt m of 
the average velocity in 1 + 1 dimensions. Data for different 
values of the force / are indicated by different symbols, rang- 
ing from 0.0204 (bottom curve) to 0.2454 (top curve). The 
system size is L — 1024 and each result was averaged over 30 
realizations of the disorder, (b) Data collapse of the veloci- 
ties according to (4.7) using the values for the DPD-3 model 
exponents (Table I). 



FIG. 7. QEW-1 model. Dependence on the tilt m of the 
average velocity in 1 + 1 dimensions. Data for different values 
of the force / are indicated by different symbols, ranging from 
0.0125 (bottom curve) to 0.075 (top curve). The system size 
is L — 1024 and each result was averaged over 30 realizations 
of the disorder. We can see clearly that A = 0. 

FIG. 8. RFIM. (a) Dependence on the tilt m of the aver- 
age velocity in 2 + 1 dimensions. Data for different values of 
the force / are indicated by different symbols, ranging from 
0.014 (bottom curve) to 0.143 (top curve). The system size 
is L 2 = 40 x 40, and A = 3 (SA regime), (b) Data collapse 
of the velocities according to (4.7) using the values for the 
RFIM exponents (Table I). 

FIG. 9. Schematic representation of the scaling behavior 
of the coefficient A of the nonlinear term for the two univer- 
sality classes. For the DPD universality class we see that A 
diverges at the depinning transition, while for the QEW uni- 
versality class it vanishes to zero at the transition. 

FIG. 10. DPD universality class, (a) Plot of the local 
width w as a function of £ for a pinned and a moving inter- 
face. The different roughness exponent is clear. The system 
size is L — 10480 and each result was averaged over 50 re- 
alizations of the disorder, (b) Plot of the local width in the 
moving phase for several values of the driving force, (c) Data 
collapse of the widths according to (5.3). (d) Data collapse 
for / = 0.0583, and for systems of different size L. Visually 
apparent is the fact that the deviations from scaling observed 
in (c) are due to finite size effects. 

FIG. 11. DPD universality class, (a) Plot of the global 
width W as a function of time the moving phase interface. 
The system size is L — 10480 and each result was averaged 
over 50 realizations of the disorder, (b) Data collapse of the 
widths according to (5.11). 

FIG. 12. Anomalous scaling of the local width for the 
QEW universality class, (a) Interfaces generated with the 
QEW-2 model, at the depinning transition, for system of sizes 
L = 128, 256, 512, 1024. We note that the local slopes increase 
with the system size, (b) Plot of the scaling of the local width 
as a function of L, for two given values of I. The fit corre- 
sponds to ctG — ap ?s 0.25. 

FIG. 13. QEW universality class, (a) Plot of the local 
width in the moving phase for several values of the reduced 
force. The system size is L = 5524 and each result was aver- 
aged over 100 realizations of the disorder, (b) Data collapse 
of the widths according to (5.3). 
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FIG. 14. QEW universality class, (a) Plot of the global 
width W as a function of time for the moving phase. The sys- 
tem size is L = 5524 and each result was averaged over 100 
realizations of the disorder, (b) Data collapse of the widths 
according to (5.11). 



FIG. 15. Simulation results for the width w(£,g) of the 
pinned interface in 1 + 1 dimensions, where g is the gradient 
in p. (a) The widths for values of the gradient ranging from 
1/70 to 1/7680. The system size is 10480 and each result is 
averaged over 30 realizations of the disorder, (b) Data col- 
lapse of these results according to the scaling form (B4), using 
the values of the exponents given in (B5). 



TABLE I. Exponents for the six models studied (see defi- 
nitions in the text). A negative value of (j> means that A — > 
when / — » [cf. Eq. (4.5)]. We argue in the text that 
the models above the horizontal line (DPD-1, DPD-2, and 
DPD-3) belong to the universality class of Eq. (3.3) and can 
be mapped, in 1 + 1 dimensions, to DP. The models below 
the line belong to the universality class of Eq. (3.1). See [56] 
for a discussion on the values of the exponents for the DPD-2 
model. 



Model 


1 + 1 dimensions 

e <t> 


2 + 1 dimensions 

6 <t> 


DPD-1 


0.64 ±0.08 0.64 ±0.08 


0.80 + 0.12 0.30 ±0.12 


DPD-2 


0.59 ±0.12 0.55 ±0.12 




DPD-3 


0.70 ±0.12 0.65 ±0.12 





QEW-1 0.26 ± 0.07 
QEW-2 0.24 ± 0.03 



FA 0.9 ±0.2 1.8 ±0.2 
RFIM SA 0.60 ±0.11 -0.70 + 0.11 

SS 0.31 + 0.08 -0.65 + 0.13 



TABLE II. Critical exponents for the two universality 
classes studied in this paper. The subscript "P" refers to 
the pinned phase, the subscript "M" denotes the quenched 
disorder regime in the moving phase, and the subscript "A" 
refers to the annealed disorder regime. The exponents z and 
za were calculated from scaling relations, while the remaining 
exponents were calculated directly in the simulations. 



Exponents 


DPD 


QEW 


Ctp 


0.63 ± 0.03 


0.92 ±0.04 


OLM 


0.75 ± 0.04 


0.92 ± 0.04 


OLA 


0.50 ±0.04 


0.46 ± 0.04 


OLG 




1.23 ±0.04 


I3p 


0.67 + 0.05 


0.85 ±0.03 


13m 


0.74 ± 0.06 


0.86 ±0.03 


Pa 


0.30 ± 0.04 


0.25 ±0.04 


z 


1.01 ±0.10 


1.45 ±0.07 


ZA 


1.67 ±0.26 


2.09 ± 0.40 


V 


1.73 ±0.04 


1.35 ±0.04 


e 


0.64 ±0.12 


0.24 ±0.03 


<fiM 


-0.12 ±0.06 


0.44 ± 0.05 


<fiA 


0.34 ± 0.06 


0.99 ±0.05 


km 


-0.11 ±0.06 


0.00 ± 0.06 


KA 


0.65. ±0.06 


1.15 ±0.06 
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